#include "../include/Multigrid_Solver.h"
#include "../include/Multigrid_Solver_1D.h"
#include "../include/Multigrid_Solver_2D.h"
#include "../include/Multigrid_Factory.h"
#include "../include/Functions_test.cpp"
#include <string>

std::string Rtype2str(int i)
{
    if (i == 1)
        return "FullWeight";
    else if (i == 2)
        return "Injection";
    else
        return "unkown";
};

std::string Itype2str(int i)
{
    if (i == 1)
        return "Linear";
    else if (i == 2)
        return "Quadratic";
    else
        return "unkown";
};

std::string Stype2str(int i)
{
    if (i == 1)
        return "WTJacobi";
    else if (i == 2)
        return "GaussSeidal";
    else
        return "unkown";
};

int main()
{
    int maxium_iteration = 50;
    std::cout << "Nonhomogeneous Possion Equation:" << std::endl;
    std::cout << "u''(P) =  (-2*sin(P[0])*sin(P[1])+sin(P[1])*cos(P[0])*cos(P[0])+sin(P[0])*cos(P[1])*cos(P[1]))*exp(sin(P[0])*sin(P[1])), and u=exp(sin(P[0])*sin(P[1])) on the boundary." << std::endl;
    int numofgrid;
    double _a[2] = {0, 0};

    VectorXd solution_temp, rhs, r_temp, r_origin, solution_exact, sol_err_box, res_err_box;
    double r_origin_norm;
    Factory2DMGSolver Solver2D;
    FactoryMGSolverBase *_p_fac_base = &Solver2D;
    Multigrid_Solver *_arena2D = _p_fac_base->CreateMGSolver();
    int flag1 = 0, flag2 = -1;
    MatrixXd _Convergence_Rate_Save;
    _Convergence_Rate_Save.setZero(8, 4);
    numofgrid = 128;
    int cnt;
    double h;
    char *Rstr, Istr;
    flag1 = 0, flag2 = -1;
    for (numofgrid = 128; numofgrid <= 1024; numofgrid *= 2)
    {
        flag1 = 0;
        flag2 += 1;
        for (int _Stype = 2; _Stype <= 2; _Stype++)
        {
            for (int _Rtype = 1; _Rtype <= 2; _Rtype++)
            {
                for (int _Itype = 1; _Itype <= 2; _Itype++)
                {
                    _arena2D->Initialize(numofgrid, Dim2_type1_df2rd, Dim2_type1_f);
                    _arena2D->setRestircType(_Rtype);
                    _arena2D->setIntpType(_Itype);
                    _arena2D->setSmoothType(_Stype);
                    solution_temp = _arena2D->get_u_origin();
                    rhs = _arena2D->get_f_origin();
                    h = _arena2D->get_h();
                    r_origin = _arena2D->ResidualSolver(solution_temp, rhs, h);
                    r_origin_norm = r_origin.lpNorm<Eigen::Infinity>();
                    solution_exact.setZero((numofgrid + 1)*(numofgrid+1));
                    for (int j = 0; j <= numofgrid; j++)
                    {
                        for (int i = 0; i <= numofgrid; i++)
                        {
                            _a[0] = 1.0 * i / numofgrid;
                            _a[1] = 1.0 * j / numofgrid;
                            solution_exact(j*(numofgrid+1)+i) = Dim2_type1_f(_a);
                        };
                    };

                    sol_err_box.setZero(maxium_iteration + 2);
                    res_err_box.setZero(maxium_iteration + 2);
                    sol_err_box(0) = solution_exact.lpNorm<Eigen::Infinity>();
                    res_err_box(0) = r_origin_norm;
                    for (cnt = 1; cnt <= maxium_iteration; cnt++)
                    {
                        solution_temp = _arena2D->Vcycle(solution_temp, rhs, h);
                        r_temp = _arena2D->ResidualSolver(solution_temp, rhs, h);
                        sol_err_box(cnt) = (solution_temp - solution_exact).lpNorm<Eigen::Infinity>();
                        res_err_box(cnt) = (r_temp).lpNorm<Eigen::Infinity>();
                        if (res_err_box(cnt) <= r_origin_norm * 1e-8)
                            break;
                    };
                    _Convergence_Rate_Save(flag1, flag2) = sol_err_box(cnt);
                    flag1++;
                    std::cout << "### 2D_VC:" << numofgrid << "," << _Stype << "," << _Rtype << "," << _Itype << "### (numofgrid,Stype,Rtype,Itype)" << std::endl;
                    std::cout << "This V-cycle with "
                              << "#n = " << numofgrid << " needs #" << cnt << " iterations for 1e-8 re_acu , when :" << std::endl;
                    printf("#Smoother               #Restricter             #Interpolater \n");
                    printf("%10s    %20s    %20s\n", Stype2str(_Stype).c_str(), Rtype2str(_Rtype).c_str(), Itype2str(_Itype).c_str());
                    std::cout << "#Detail:" << std::endl;
                    printf("%17s %15s %15s %15s %15s\n", ((std::string) "Relative_Accuracy").c_str(), ((std::string) "Residual").c_str(), ((std::string) " Ratio_Residual").c_str(), ((std::string) "SolutionError").c_str(), ((std::string) "Ratio_error").c_str());
                    for (int i = 1; i <= cnt; i++)
                    {
                        printf("%2d %14.2e %14.2e %14.2e %14.2e %14.2e\n", i, res_err_box(i) / res_err_box(0), res_err_box(i), res_err_box(i) / res_err_box(i - 1), sol_err_box(i), sol_err_box(i) / sol_err_box(i - 1));
                    };
                    std::cout << "Wait for about 20 seconds ..." << std::endl;
                };
            };
        };
    };
std::cout<<"EOF"<<std::endl;
    // MatrixXd _Generater_convergence;
    // _Generater_convergence.setZero(8, 3);
    // for (int j = 0; j <= 2; j++)
    // {
    //     for (int i = 0; i <= 7; i++)
    //     {
    //         _Generater_convergence(i, j) = -log2(abs(_Convergence_Rate_Save(i, j + 1) / _Convergence_Rate_Save(i, j)));
    //     };
    // };

    // std::cout << "##Convergence Level Table" << std::endl;
    // MatrixXd A(1, 4);
    // A << 0, 128, 256, 512;
    // MatrixXd B(5, 1);
    // B << 0, 112, 122, 212,222;
    // MatrixXd C(5, 4);
    // C.row(0) = A;
    // C.col(0) = B;
    // for (int j = 1; j <= 3; j++)
    //     for (int i = 1; i <= 4; i++)
    //         C(i, j) = _Generater_convergence(i - 1, j - 1);
    // std::cout << C << std::endl;
    // std::cout << "PS:The first col is the combination code of IntpType,RestrictType and SmootheType." << std::endl;
    // std::cout << "Interpolation Function: Linear 　= 　1；Quadratic 　= 　2；" << std::endl;
    // std::cout << "Restriction Function:   FullWeighting 　= 　1； Injection 　= 　2；" << std::endl;
    // std::cout << "Smoothing Function: WeightedJacobi ,omega = 2/3　= 　1;   GaussSeidel 　= 　2；" << std::endl;
};